#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>
#include <fstream>
#include "Spline.hpp"
#include "Function.hpp"

class F1 : public Function {
public:
    double operator() (double x) const override {
        return 1 / (1 + 25 * x * x); // Example function
    }
};

double maxNormError(const std::vector<double>& errors) {
    double maxError = 0.0;
    for (double error : errors) {
        if (std::abs(error) > maxError) {
            maxError = std::abs(error);
        }
    }
    return maxError;
}

void computeErrorsAndConvergenceRates(const std::vector<int>& Ns) {
    F1 func;
    std::vector<double> linearErrors;
    std::vector<double> cubicErrors;
    std::ofstream linearOutfile("output/A_Linear.txt");
    std::ofstream cubicOutfile("output/A_Cubic.txt");
    for (int N : Ns) {
        std::vector<double> knots(N);
        std::vector<double> values(N);
        double h = 2.0 / (N - 1);
        for (int i = 0; i < N; ++i) {
            knots[i] = -1 + i * h;
        }

        // Generate linear PPForm and CompleteCubicPPForm splines
        LinearPPForm linearSpline(knots, func);
        CompleteCubicPPForm cubicSpline(knots, func);

        // Calculate errors
        std::vector<double> linearSplineErrors(N - 1);
        std::vector<double> cubicSplineErrors(N - 1);
        for (int i = 0; i < N - 1; ++i) {
            double midPoint = (knots[i] + knots[i + 1]) / 2.0;
            double exactValue = func(midPoint);
            double linearSplineValue = linearSpline.evaluate(midPoint);
            double cubicSplineValue = cubicSpline.evaluate(midPoint);
            linearSplineErrors[i] = exactValue - linearSplineValue;
            cubicSplineErrors[i] = exactValue - cubicSplineValue;
        }

        // Plot functions
        for (double x = knots.front(); x <= knots.back(); x += 0.01) {
            double linearY = linearSpline.evaluate(x);
            double cubicY = cubicSpline.evaluate(x);
            linearOutfile << x << " " << linearY << std::endl;
            cubicOutfile << x << " " << cubicY << std::endl;
        }
        linearOutfile << std::endl;
        cubicOutfile << std::endl;

        // Record maximum errors
        double maxLinearError = maxNormError(linearSplineErrors);
        double maxCubicError = maxNormError(cubicSplineErrors);
        linearErrors.push_back(maxLinearError);
        cubicErrors.push_back(maxCubicError);
    }

    // Output errors to file
    std::ofstream errorOutfile("output/Aerrors.txt");
    for (int i = 0; i < Ns.size(); ++i) {
        errorOutfile << std::log(Ns[i]) << " " << std::log(linearErrors[i]) << std::endl;
    }
    errorOutfile << std::endl;
    for (int i = 0; i < Ns.size(); ++i) {
        errorOutfile << std::log(Ns[i]) << " " << std::log(cubicErrors[i]) << std::endl;
    }
    linearOutfile.close();
    cubicOutfile.close();
    errorOutfile.close();
}

int main() {
    std::vector<int> Ns = {6, 11, 21, 41, 81};
    computeErrorsAndConvergenceRates(Ns);
    return 0;
}
